Packages / R-options:

.libPaths('~/R/x86_64-pc-linux-gnu-library/4.1/')
suppressPackageStartupMessages({
  library(plotly)
  library(gaston)
  library(dplyr)
})
# load engine's functions
sapply(FUN = source,
       X = list.files("../../src",
                      pattern = ".R$",
                      full.names = T))

#  R options
options(max.print = 20)

1 Get data

1.1 Inputs

genoFile <- '../../data/geno/breedGame_phasedGeno.vcf.gz'
crossTableFile <- '../../data/crossingTable/breedGame_small_crossTable.csv'
SNPcoordFile <- '../../data/SNPcoordinates/breedingGame_SNPcoord.csv'
markerEffectsFile <- '../../data/markerEffects/breedGame_markerEffects.csv'

1.2 Predictions

predictions <- calc_progenyBlupEstimation(
  genoFile = genoFile,
  crossTableFile = crossTableFile,
  SNPcoordFile = SNPcoordFile,
  markerEffectsFile = markerEffectsFile,
  # outFile = outFile
)
predictions$Cross <- paste0(predictions$ind1, '_X_', predictions$ind2)

1.3 Crossing simulation

set.seed(1234)
markerEffects <- readMarkerEffects(markerEffectsFile)
crossSimOutFile <- crossingSimulation(
  genoFile = genoFile,
  crossTableFile = crossTableFile,
  SNPcoordFile = SNPcoordFile,
  nCross = 100)
simulation <- gaston::read.vcf(crossSimOutFile)

## Warning in gaston::read.vcf(crossSimOutFile): NAs introduced by coercion

## Warning in gaston::read.vcf(crossSimOutFile): Some unknown chromosomes id's (try
## to set convert.chr = FALSE)

simulation <- gaston::as.matrix(simulation)
simulation <- simulation[, row.names(markerEffects$SNPeffects)]

simulation <- data.frame(
  simBV = as.vector(simulation %*% as.matrix(markerEffects$SNPeffects)),
  Cross = as.factor(sub('-\\d+$', '', row.names(simulation)))
)

2 Simulation VS Predictions

p <- plot_ly() 

# add boxplot of simulated BV
p <- p %>% add_boxplot(data = simulation,
                       x = ~Cross,
                       y = ~simBV,
                       name = "Boxplot Simulated BV")
# add jittered markers of simulated BV
p <- p %>% add_boxplot(data = simulation,
                       x = ~Cross,
                       y = ~simBV,
                       line = list(color = 'rgba(0,0,0,0)'),
                       marker = list(color = 'rgb(31, 119, 180)'),
                       fillcolor = 'rgba(0,0,0,0)',
                       name = "Markers Simulated BV",
                       boxpoints = "all",
                       pointpos = 0,
                       jitter = 1,
                       hoverinfo = 'text',
                       text = apply(simulation, 1, function(l) {
                         paste(names(l), ":", l, collapse = "\n")
                       }))


# add predicted BV
p <- p %>% add_markers(
  data = predictions,
  x = ~Cross,
  y = ~blup_exp,
  marker = list(color = 'rgb(255,0,0)'),
  name = "Predicted BV",
  error_y = ~ list(array = sqrt(blup_var),
                   type = 'data',
                   color = '#000000'),
  hoverinfo = 'text',
  text = apply(predictions, 1, function(l) {
    paste(names(l), ":", l, collapse = "\n")
  })
)

p

2.1 Test observed mean

alpha <- 0.05
beta <- 0.95
pred_vs_sim <- full_join(simulation, predictions, by = "Cross") %>% 
  filter(Cross != "Coll0659_X_Coll0425") %>% 
  group_by(Cross) %>%
  summarise(obs_mean =  mean(simBV),
            blup_exp = unique(blup_exp),
            p.value = t.test(x = simBV, mu = unique(blup_exp))$p.value,
            delta = power.t.test(n = 100, power = beta, sd = sqrt(unique(blup_var)), sig.level = alpha)$delta)
pred_vs_sim$adj.p.val <- p.adjust(pred_vs_sim$p.value, method = "bonferroni")
pred_vs_sim

idLine <- data.frame(mean = c(min(c(pred_vs_sim$blup_exp, pred_vs_sim$obs_mean)) - 1,
                                max(c(pred_vs_sim$blup_exp, pred_vs_sim$obs_mean)) + 1),
                     var = c(min(c(pred_vs_sim$blup_var, pred_vs_sim$sim_var)) - 1,
                                max(c(pred_vs_sim$blup_var, pred_vs_sim$sim_var)) + 1))

## Warning: Unknown or uninitialised column: `blup_var`.

## Warning: Unknown or uninitialised column: `sim_var`.

## Warning in min(c(pred_vs_sim$blup_var, pred_vs_sim$sim_var)): no non-missing
## arguments to min; returning Inf

## Warning: Unknown or uninitialised column: `blup_var`.

## Warning: Unknown or uninitialised column: `sim_var`.

## Warning in max(c(pred_vs_sim$blup_var, pred_vs_sim$sim_var)): no non-missing
## arguments to max; returning -Inf

rmse <- sqrt(mean((pred_vs_sim$blup_exp - pred_vs_sim$obs_mean)^2))

RMSE = 0.2700433

#plot mean
plot_ly(type = "scatter",
        mode = "markers",
        data = pred_vs_sim,
        x = ~blup_exp,
        y = ~obs_mean,
        name = "Observed mean vs prediction",
        hoverinfo = 'text',
        text = apply(pred_vs_sim, 1, function(l) {
          paste(names(l), ":", l, collapse = "\n")
        })
) %>%
  add_lines(inherit = FALSE,
            x = idLine$mean,
            y = idLine$mean,
            name = "Identity line")

2.2 Test normality

pred_vs_sim <- full_join(simulation, predictions, by = "Cross") %>% 
  filter(Cross != "Coll0659_X_Coll0425") %>% 
  group_by(Cross) %>%
  summarise(obs_var = var(simBV),
            blup_var = unique(blup_var),
            p.value = ks.test(simBV, "pnorm", mean = unique(blup_exp), sd=sqrt(unique(blup_var)))$p.value)
pred_vs_sim$adj.p.val <- p.adjust(pred_vs_sim$p.value, method = "bonferroni")
pred_vs_sim

Appendix

Session Information (click to expand)

## Document generated in:
## Time difference of 26.87301 secs
## 
## CPU: AMD Ryzen 5 3600X 6-Core Processor
## Memory total size: 32.79254 GB
## 
## 
## Session information:
## R version 4.1.1 (2021-08-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Pop!_OS 22.04 LTS
## 
## Matrix products: default
## BLAS:   /usr/local/lib/R/lib/libRblas.so
## LAPACK: /usr/local/lib/R/lib/libRlapack.so
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] dplyr_1.0.7        gaston_1.5.7       RcppParallel_5.1.5 Rcpp_1.0.8        
## [5] plotly_4.10.0      ggplot2_3.3.5     
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.1.1    xfun_0.29           bslib_0.3.1        
##  [4] memuse_4.2-1        purrr_0.3.4         splines_4.1.1      
##  [7] lattice_0.20-44     colorspace_2.0-2    vctrs_0.3.8        
## [10] generics_0.1.2      htmltools_0.5.2     viridisLite_0.4.0  
## [13] vcfR_1.12.0         mgcv_1.8-36         yaml_2.2.2         
## [16] utf8_1.2.2          rlang_1.0.1         jquerylib_0.1.4    
## [19] pillar_1.7.0        glue_1.6.1          withr_2.4.3        
## [22] DBI_1.1.2           lifecycle_1.0.1     stringr_1.4.0      
## [25] munsell_0.5.0       gtable_0.3.0        htmlwidgets_1.5.4  
## [28] breedSimulatR_0.3.2 evaluate_0.14       knitr_1.37         
## [31] permute_0.9-7       fastmap_1.1.0       crosstalk_1.2.0    
## [34] parallel_4.1.1      fansi_1.0.2         pinfsc50_1.2.0     
## [37] scales_1.1.1        vegan_2.6-2         jsonlite_1.7.3     
## [40] digest_0.6.29       stringi_1.7.6       grid_4.1.1         
## [43] cli_3.1.1           tools_4.1.1         magrittr_2.0.2     
## [46] sass_0.4.0          lazyeval_0.2.2      tibble_3.1.6       
## [49] cluster_2.1.2       ape_5.6-2           crayon_1.4.2       
## [52] tidyr_1.2.0         pkgconfig_2.0.3     Matrix_1.3-4       
## [55] MASS_7.3-55         ellipsis_0.3.2      data.table_1.14.2  
## [58] assertthat_0.2.1    rmarkdown_2.11      httr_1.4.2         
## [61] rstudioapi_0.13     R6_2.5.1            nlme_3.1-152       
## [64] compiler_4.1.1

shiny::HTML('&nbsp;
<!-- Add icon/font library -->
<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.13.0/css/all.min.css">
<link rel="stylesheet" href="https://fonts.googleapis.com/css?family=Lato">

<div class="footer" style="font-family: Lato">
    <hr />
    <p style="text-align: center;"><a href="https://github.com/juliendiot42">Julien Diot</a></p>
    <p style="text-align: center;"><span style="color: #808080;"><em>juliendiot@ut-biomet.org</em></span></p>

<!-- Add font awesome icons -->
    <p style="text-align: center;">
        <a href="https://github.com/juliendiot42" class="fab fa-github"></a>
        <a href="https://www.linkedin.com/in/julien-diot-949592107/" class="fab fa-linkedin"></a>
        <a href="https://orcid.org/0000-0002-8738-2034" class="fab fa-orcid"></a>
        <a href="https://keybase.io/juliendiot" class="fab fa-keybase"></a>
    </p>
</div>
&nbsp;')

Crossing simulation vs Predictions

Julien Diot

2022/10/04

Packages / R-options:

1 Get data

1.1 Inputs

1.2 Predictions

1.3 Crossing simulation

2 Simulation VS Predictions

2.1 Test observed mean

2.2 Test normality

Appendix